/**
 * @file TimeIntegrator.h
 * @brief 
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-04-26
 * 
 * @copyright Copyright (c) 2021  Linhuo
 * 
 * @par 修改日志:
 * <table>
 * <tr><th>Date       <th>Version <th>Author  <th>Description
 * <tr><td>2021-04-26 <td>1.0     <td>Linhuo     <td>内容
 * </table>
 */
#ifndef _lin_TimeIntegrator_
#define _lin_TimeIntegrator_

#include <fstream>
#include <string>
#include <cmath>
#include <cstdlib>

#include <iostream>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Dense>

#include "func_detail.h"

class TimeIntegrator
{
protected:
    int s;
    int p;
    double dt;
    std::string name; 
    Eigen::MatrixXd u;/// saving temp data of u
    Eigen::MatrixXd fu;/// saving temp data of fu
    Eigen::VectorXd rough_sol;///saving rough solution

public:
    TimeIntegrator(){};

    void set_dt(double _dt){
        dt = _dt;
    };

    double get_dt() const
    {
        return dt;
    };

    void list();

    virtual void load_scheme(std::string _name)=0;

    void initial_value(Eigen::VectorXd _u);    

    virtual void One_step_forward()=0;

/**
 * @brief Served as the engine of ODEsolver.
 * @param  T                The given period.
 */
    void solution_generator(double T);

/**
 * @brief Served as the shell of ODEsolver. Return a numerical solution, along with generating a Matlab code
 * that is easy to plot.
 * @param  _name            Method name along with accuacy,such as "AB2" "ERK4".
 * @param  _u               Initial value.
 * @param  _T               The given period.
 * @return Eigen::VectorXd 
 */
    Eigen::VectorXd Create_as_uwish(std::string _name, Eigen::VectorXd _u,double _T);

    Eigen::VectorXd Create_as_uwish(std::string _name, Eigen::VectorXd _u,double _T,double _dt);
};

Eigen::VectorXd TimeIntegrator::Create_as_uwish(std::string _name, Eigen::VectorXd _u,double _T)
{
    load_scheme(_name);
    initial_value(_u);
    solution_generator(_T);
    return rough_sol;
};

Eigen::VectorXd TimeIntegrator::Create_as_uwish(std::string _name, Eigen::VectorXd _u,double _T,double _dt)
{
    load_scheme(_name);
    initial_value(_u);
    set_dt(_dt);
    solution_generator(_T);
    return rough_sol;
};

void TimeIntegrator::solution_generator(double T)
{   
    double t = 0.0;
    std::string buf="Mat_output_"+name+".m";
    char file_name[20];
    strcpy(file_name,buf.c_str());
    std::ofstream fout(file_name);
    fout << "u_"<< name << "=[" << std::endl;
    while (t < T)
    {
	One_step_forward();
	for (int i = 0; i < 5; i++)
	    fout << u(i,s-1) << ", ";
    fout << u(5,s-1) << " ;" << std::endl;
	t += dt;
    }
    fout << "];" << std::endl;
    fout <<"plot(u_"<< name << "(:,1),u_" << name << "(:,2))" <<std::endl;
    fout.close();
    rough_sol.resize(6);
    rough_sol = u.col(s-1);
}

void TimeIntegrator::initial_value(Eigen::VectorXd _u0)
{   u.resize(6,s);
    fu.resize(6,s);
    u.col(0) = _u0;
    fu.col(0) = (func)(_u0);
    Eigen::VectorXd y_1,y_2,y_3,y_4;

    for(int i=0;i < s - 1; i++)
    {
    // std::cout << "..(initializing)" << std::endl;
        y_1 = (func)(u.col(i));
        y_2 = (func)(u.col(i) + 0.5 * dt * y_1 );
        y_3 = (func)(u.col(i) + 0.5 * dt * y_2 );
        y_4 = (func)(u.col(i) + dt * y_3);
        u.col(i+1) = u.col(i) + (dt / 6.0) * (y_1 + 2 * y_2 + 2 * y_3 + y_4 );
        fu.col(i+1) = (func)(u.col(i+1));
    };
};

void TimeIntegrator::list()
{
    std::cout << "s = " << s << std::endl;
    std::cout << "p = " << p << std::endl;
    std::cout << "dt = " << dt << std::endl;
    std::cout << "name = " << name << std::endl;
};

#else
//DO NOTHING
#endif
